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For weakly non ergodic systems, the probability density function of a time average observable O 



is f a (o) = -ilim E _ Im- 



where Oi is the value of the observable when the 



system is in state i = !,■■■ L. p° q is the probability that a member of an ensemble of systems occupies 
state i in equilibrium. For a particle undergoing a fractional diffusion process in a binding force 
field, with thermal detailed balance conditions, p^ q is Boltzmann's canonical probability. Within 
the unbiased sub-diffusive continuous time random walk model, the exponent < a < 1 is the 
anomalous diffusion exponent (x 2 ) ~ t a found for free boundary conditions. When a — * 1 ergodic 
statistical mechanics is recovered lim a _,i f a ((D) =5(0-(0)). We briefly discuss possible physical 
applications in single particle experiments. 

PACS numbers: 05.70.Ln, 05.20.Gg, 05.40.Fb 



I. INTRODUCTION 

An ensemble of non-interacting one dimensional Brow- 
nian particles in the presence of a binding potential field 
V(x) reach a thermal equilibrium described by Boltz- 
mann's canonical law p eq (x) — exp[— V(x)/T]/Z where 
T is the temperature (units fc& = 1) and Z is the nor- 
malizing partition function. With this law we may calcu- 
late ensemble averages for example (x) = f_ xp eq (x)dx. 
On the other hand from the trajectory of a single particle 
x(t) we may construct the time average x = Jl* x(t')dt'/t. 
For ergodic motion the time and ensemble averages are 
identical in the limit of long measurement times t — > oo. 
What is the Physical meaning of a long measurement 
time? Brownian dynamics in a finite interval, is char- 
acterized by a finite relaxation time which is the time 
scale on which particles reach thermal equilibrium. For 
the simplest case of Brownian motion between two re- 
flecting walls, with a system of size I, dimensional anal- 
ysis gives a relaxation time of the order l 2 /D where D 
is the diffusion constant. A second time scale, the aver- 
age time between jump events (r) is more microscopical. 
The latter is related to D with the Einstein relation [l| 
D = ((Ax) 2 ) /2(r) where ((Ax) 2 } is the variance of jump 
lengths. Ergodicity is almost trivial when the measure- 
ment time is much longer than these two time scales. For 
example for a Brownian motion in a Harmonic field. 

On the other hand anomalous diffusion and transport, 
is characterized in many cases by a diverging relaxation 
time and a diverging microscopical time scale (r) . For ex- 
ample unbiased sub-diffusion is characterized by a mean 
square displacement (a; 2 ) oc t a and < a < 1. The 
reader immediately realizes that the diffusion constant 
is zero, in the sense that \im t ^ (x >(x 2 ) /t = 0, hence the 
mentioned relaxation time I 2 / D is infinite even when the 
system size / is finite. Indeed according to the contin- 
uous time random walk (CTRW) model @, i, H i, i], 
anomalous sub-diffusion is found when waiting times be- 
tween jumps diverge (r) — * oo, which is related to the 
Scher-Montroll power law waiting time probability den- 



sity function (PDF) iP(t) oc t~( 1+q ) 2]. For such scale 
free anomalous diffusion the relaxation time and the av- 
eraged sojourn time (r) are infinite, and ergodicity is 
broken weakly. 

Strong ergodicity breaking is found when a system 
is divided into inaccessible regions of its phase space. 
Namely a particle or a system starting in one region can- 
not explore all other regions due to some non-passable 
barrier (e.g. in the micro-canonical case, by regions we 
mean sections on the constant energy surface) . Bouchaud 
[3] introduced the profound concept of weak ergodicity 
breaking in the context of glass dynamics, which in turn 
is related to infinite ergodic theory Q investigated by 
Mathematicians using a dynamical approach. In weak 
ergodicity breaking the phase space is not broken into in- 
accessible regions. Instead due to the power law sticking 
times the dynamics is non-stationary and non-ergodic. 
Since the ergodic hypothesis is the pillar on which sta- 
tistical mechanics is built, but at the same time also 
long tailed power law distributions of trapping times are 
very common in the description of Physical behaviors 
0- HI El IE @| , it is natural to investigate the non-ergodic 
properties of systems whose stochastic dynamics are gov- 
erned by such anomalous statistics. Previously weak er- 
godicity breaking was investigated for: blinking quan- 
tum dots 0, [HI EE E3] , Levy walks_|3 occupation time 
statistics of the CTRW model Il5j. fractional Fokker- 
Planck equation [l6| , deterministic one dimensional maps 
17, 18, 19], numerical simulations of fractional transport 
in a washboard potential |2 C| and in vivo gene regulation 
by DNA-binding proteins [2l|. Recently a relation be- 
tween statistics of weak ergodicity breaking and statistics 
of non-self averaging in models of quenched disorder was 
found [HJ . Hence it is timely to present a general statis- 
tical mechanical framework for weak-ergodicity breaking. 

In this manuscript we investigate the distribution of 
time averaged observables for weak ergodicity breaking. 
We explore the relations between ensemble averages and 
fluctuations of time averages. And investigate the tran- 
sition from the localization limit a — > to the usual 
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ergodic behavior found for a — > 1. Specific examples 
for the distribution of x for a particle undergoing a sub- 
diffusive CTRW in a potential V(x) are worked out is 
detail. In the second part of the paper we derive our 
main results using a generalized CTRW approach. Wc 
investigate models with a single waiting time PDF and 
more general models with several types of such PDFs. A 
brief summary of some of our results was recently pub- 
lished in 

The theory of weak-ergodicity breaking is mathemati- 
cally related to the arcsine distribution [24]) [2^, [26| , and 
to its extensions [13, . Consider a normal Brownian 
motion with free boundary conditions in one dimension 
x(t) = rj(t) where rj(t) is Gaussian white noise. The time 
t + spent by the particle in x > is called an occupation 
time. Naive expectation is that the single particle will 
occupy x > for half of the measurement time t, when 
the latter is long. Instead the PDF of the occupation 
fraction is 



/ 



t+ 



t J n^(t+/t)(l-t+/t) 



(1) 



This arcsine law is related to the well known PDF of 
first-passage times from x to the origin, for simple Brow- 
nian motion. The latter PDF decays i -3 / 2 for long first- 
passage times [25j |. and the averaged return time is in- 
finite. Roughly speaking, during the dynamics of the 
particle, it will usually occupy either the domain x > 
or the domain x < for a duration which is of the order 
of the measurement time. Thus the arcsine PDF Eq. ([1]) 
has a U shape. This behavior is found because the Brow- 
nian motion was assumed to be unbounded. If we add 
reflecting walls on x = 1/2 and x = —1/2 the dynamics 
will be ergodic and in the long time limit the particle 
spends half of the time in (0, 1/2) and the other half in 
(—1/2,0). The theory of weak ergodicity breaking, pre- 
sented in this manuscript, is mathematically related to 
the arcsine law and is based on Levy statistics. 

This paper is organized as follows. In Sec. |TT] distribu- 
tion of time averaged observables for weakly non-ergodic 
systems is presented using general arguments not spe- 
cific to a model. Properties of this distribution are in- 
vestigated in Sec. UTT] and in Sec. |IV] the example of 
the distribution of x for a particle undergoing fractional 
dynamics in a binding potential is worked out in detail. 
In Sec. [V]we derive our main result using a CTRW ap- 
proach, thus further justifying assumptions made in Sec. 
ITTl Numerical simulations of x obtained from the CTRW 



process are compared with analytical theory in sub-Sec. 
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II. DISTRIBUTION OF TIME AVERAGED 
OBSERVABLES 

We consider a system with L states and label them 
with an index i = 1, ...L. A time average of a Physical 
observable O is made. If the system is in state i the 
Physical observable attains the value Oi. The time the 
system spends in state i is ti and is called a residence 
time or an occupation time. The time average of the 
Physical observable is 



q _ Ei=l t&i 



Ei=l *i 



(2) 



and min{Oi} < O < max{0;}. As mentioned in the 
introduction many physical systems, in their stochastic or 
deterministic dynamics, are known to be characterized by 
power law sojourn times in the states of the system 0, H, 
0) S @ • We assume that the occupation time tj is a sum 
of many such sojourn times. If the state i is visited many 
times, and the sojourn times are independent identically 
distributed random variables, Levy's limit theorem will 
describe the statistics of the residence times fj in the 
limit of long measurement time. Hence we argue that 
the PDF of ti is a one sided Levy PDF l at p^(ti) whose 
Laplace transform is 



L, P ^ (U) ex.p(-uU)dti = exp (-p cq u a ) (3) 



and < a < 1 [29(. Later we find that p cq is the proba- 
bility that a member of an ensemble of systems occupies 
state i in equilibrium. For the CTRW with thermal de- 
tailed balance conditions p° q is Boltzmann's probability 
of finding the system in state i, p^ q = ex-p(—Ei/T)/Z as 
we will show later. 

The following generating function [26| is a useful tool 



/a(0 



1 



^(-i) fe <oV 



(4) 



fe=0 



Our main aim is to find the PDF of the time averaged 
observable 



f a (0)=(5 fe-IgL^Y 



■— lim Im( 

7T e^O 



o 



lim Im= — 

7T e^O O + 



1 - 



O+ie 



(•5) 
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using Eq. Q 



f a (0) =-ilimIm=^— /„ 



7T e-*o £> + ie \ O + ie 



1 



We now find the generating function Eq. ^ and invert it using Eq. ^ to obtain f a (O) . 
The generating function is rewritten 

( y L ^ 



(0) 



00 />00 poo 

ds dt I dtil a p oq (ti) 
Jo Jo 



Using a well known presentation for the delta function in Eq. ([7]) 



(7) 



2tt 



dke 



changing variables kt — k and using Eq. ([3]) we find 

/«(£) = / ^ / !P / dsexp 
JO J -00 27r Jo 

We again change variables k = k/t and s = s/t and obtain 



ijb _ s _ y >— 



(8) 



(9) 



00 poo 



dk 



/*(£)=/ dt / / dStexp 

/O J-oo 27r Jo 



jfct - it - ^p° q (ik + OiCs) a 



i=l 



This equation is rewritten using a simple trick 



/« (0 = / dt 
Jo 



00 />oo 



dk f°° , I d 
— / ds < — — exp 

00 2irj ds 



(10) 



a V (tfc + O^S) 0-1 0^ exp 



(^-S)t-2jp° q (i£: + 0^S) c 



(11) 



Integration over t gives a simple pole l/(ik — s), using 
Cauchy integral formula to integrate over fc, and then 
solving two trivial integrals yields 



fa (0 



EtiP- q (i + ^) Q 



Inverting Eq. 



12| using Eq. (J6]) we find 



7T e^0 



. \ t 



(12) 



(13) 



EtiPr(S-Oi + ie)' 
This is a very general formula for the distribution of 



time averaged observables for weakly non-ergodic sys- 
tems. As we show later, within the CTRW model, a 
is the anomalous diffusion exponent. For < a < 1 
we use lim e ^o (O - O, + ie) a = {O - J | t V^ Q where 
7T if ~0< Oi 

if O > Oi 



and Eq. (|T3|) becomes 



f a (o) 
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sin Tra I<_ x (0) l| (0) + (O) I< (O) 



/| (O) + [J< (O)] 2 + 2/| (O) I< (O) cos , 



with 



and 



o<o, 



0>Oi 



(14) 
(15) 

(16) 



Notice the L divergences of f a (0) when = 0^ due to 
the J ( ^_ 1 (0) term in the numerator of Eq. I|14p. This 
behavior is caused by long sticking times in a state of 
the system, on a time scale which is of the order of the 
measurement time. 



III. STATISTICS OF WEAK ERGODICITY 
BREAKING 

In this Sec. we investigate properties of Eqs. (|12|13|) . 

A. Limits a — > and a — * 1 

In the limit a — > 1 we recover usual ergodic behavior. 
From Eq. (|lc 



(17) 



Using the normalization condition J2i=iPi = 1 an< ^ the 
ensemble average 



i=l 



we have ergodic behavior 

A(O) =5(0-(0)) 



(18) 



(19) 



in the sense that the ensemble average is equal to the 
time average. Note that already Eq. in the limit 
a — ► 1 indicates ergodicity since the residence time is not 
a fluctuating quantity (i.e. the one sided Levy PDF is 
a delta function when a = 1). In the opposite limit of 
a — > we find 

_ 1 L _ 

lim / Q (0) = lim Im V (5 - C 4 + ie) 1 (20) 

<=1 

hence 



This describes a localization behavior where the system 
is stuck in one of the states for the whole duration of the 
observation time, which is the expected behavior when 
a -> 0. 



B. Lamperti Statistics of the Occupation fraction 

Let the Physical observable be Oi — 1 when i = 1, I 
where I < L, otherwise Oi — 0. Hence the time average 
in this case is the occupation fraction 



(22) 



which is the fraction of time spent by the system in the 
observation domain i — 1, • • • , I. Clearly < O < 1 in 
this case. Using Eq. (|13[) a straight forward calculation 
gives 



1 



/- (?) = - 



K[(l-0)0] a 1 sin Tra 



t ft 2 (1 - O) + O + 2TZ [(1 - O) O] a cos na 

(23) 

This is the Lamperti PDF [2JJ which is a natural general- 
ization of the arcsine distribution [the case a = 1/2, 1Z — 
1 Eq. (TJJ)]. The PDF Eq. (23]) has found several ap- 



plications in non-ergodic s yste ms mentioned in the in- 
troduction pH Q [3 EE lH, HI and recently for the 
non-self averaging properties of the quenched trap model 
(22l ] . The parameter 1Z in Eq. ([23|) is called the asymme- 
try parameter and is given by 



K = 



eq 



cq 



(24) 



where p c ^j = J2i=i Pi* ^ s the probability in ensemble 
sense to be in the observation domain. Since the ob- 
servable attains two values Oi = 1 or Oi = the PDF 
Eq. (|23p has two divergences on O = 1 and = which 
is a special case of the more general rule discussed after 

Eq. (nnj). 



C. Low order moments of time averaged 
observables 



From the moment generating function / Q (£) we can 
obtain moments of the time averages O 



(^> = (-i)"^/-(Olc=o- 



(25) 



Using Eq. (JT2J) we find 



lim f a (p) ^pfSiO-Oi) 



(21) 



(26) 
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The average (....) is over an ensemble of realizations. 
If the ensemble reaches an equilibrium then obviously 
(O) = (O) which is time independent. Hence the p^ q s in 
Eq. (f2"6"| are the probabilities that a member of an ensem- 
ble of systems occupies state i when the ensemble reaches 
an equilibrium. This justifies our original assumption 
that the p° q s in Eq. §5§ are population fractions. Namely 
the p° q s can be in principle measured by letting many in- 
dependent systems (or many non-interacting particles) 
evolve and then in the long time limit p° q is the number 
of systems in state i over the total number of systems 
when the latter is large. Using Eqs. (|12I25[) the fluctua- 
tions are given by 

(O 2 ) (Of = (1 - a) ((O 2 ) - (O) 2 ) (27) 



where (O 2 ) = T,f=iPT(°i) 2 - E( l- <E3 S ives a sim P le 
relation between fluctuations of time averages and fluc- 
tuations of ensemble averages. Once again when a = 1 
the fluctuations of the time average Eq. (|27p vanish, in- 
dicating ergodic behavior. Relations between cumulants 
of time average observables and cumulants of ensemble 
averages are found in a similar way. For the third cumu- 
lant 



C 3 (0) = (O 3 ) - 3(0 2 )(6) + 2{6) 3 = ^(2 - a)(l - a) ((O 3 ) - 3{0 2 )(0) + 2(0) 3 ) = 1(2 - a)(l - a)C 3 (0) (28) 



and the fourth cumulant 



C 4 (0) = (O 4 ) - 4(0 3 )(0) - 3(0 2 ) 2 + 12{0 2 ){0) 2 - 6(0)- 



CI -ci [' (3 01 f a) ((O 4 ) - 4(0 3 )(0)) + (6 - 6a + a 2 ) ((O) 2 (2(0 2 ) - (O) 2 )) - 1(6 - a)(l - a){0 2 ) 2 



r 



i29) 



Low order moments of time averaged observables can be 
expressed using the cumulants Eqs. (|27l28l29p the re- 
sults remain as cumbersome as the expressions in Eqs. 
(|28l29p . When the odd moments of the ensemble average 
are equal zero (0 2n ^ 1 ) = n = 1,2, • • corresponding 
to examples we investigate in Sec. IIV1 the expressions 
for moments are simpler. Odd moments of the time av- 
erage observable are equal zero, as expected. The second 
moment is (O ) = (1 — a)(0 2 ) and the fourth moment 



r > = { - 



a)(2-a)(l-a) 



(O 4 



a(l - a) 2 



(O 2 ) 2 . 
(30) 



D. Correlations Between Occupation Fractions 

(PiPk) 

The correlations between the occupation fractions p l = 
ti/t and p k = tk/t where t = $^ i=1 it are now briefly 
investigated. We use the L dimensional vector £ = 
£2, •••£.£,} and the L dimensional generating function 

S ° ® = ( TTEb5->- (31) 

Related multidimensional arcsine distribution of occupa- 
tion fractions were investigated in (28|. Using the sub- 



stitution £0j — > £, in Eq. (g]) it is easy to see using Eq. 



fjc 



(0 



a-l 



This multidimensional generating function yields 



■9a (I) ||- = 



ft 



(32) 



(33) 



and similarly by taking the second order derivative of 
9a (0 w ith respect to & 



(^>-fo> 2 = (1-^(1-^) 



(34) 



Identical results can be obtained using the Lamperti PDF 
Eq. (j2"3")l . More interesting is to notice the correlations 
between occupation fractions for example 



/— \- 1 9 9 ~ 
{PlPk) ~ 2dfkWi 9 ° 



(I) l^o 



for I ^ k. Using Eq. (32 



(35) 



(36) 



We see that when a — * 1 the occupation fractions are 
uncorrelated since (piPk) — (Pi) (p^) — and one can show 
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that they are independent random variables. When a — > 

the system occupies one state for practically the whole 
duration of the measurement hence either p l ~ 1 and 
then obviously p k ~ or the opposite situation is found, 
or both occupation fractions are zero (if L > 2). In any 
case clearly the product p{p k is zero when a — > and 

1 k as we have indeed found in Eq. ([551) . 



where 



and 



4< (a?) 



4> (a?) 



da;p eq (x)|a;-x| a 



f dxp CCL {x)\x - x\ a 

J — OO 



(41) 



(42) 



IV. DISTRIBUTION OF x 

We now consider a particle undergoing stochastic frac- 
tional dynamics in a binding field. The fractional 
Fokker-Planck equation HH describes anomalous 
sub-diffusion and relaxation close to thermal equilibrium 
using fractional calculus 



d a p(x,t) 



D n 



d 2 d F(x) 
dx 2 ~ dx T 



p{x,t) 



(37) 



where D a is the fractional diffusion coefficient, and 
F(x) — —dV(x)/dx is the force. Eq. [)37p reduces to 
the usual Fokker-Planck equation when a = 1. The 
fractional Fokker-Planck equation was derived from the 
sub-diffusive continuous time random walk [33| which 
is the stochastic process we have in mind. In the ab- 
sence of the force field and for free boundary conditions 
(x 2 ) = 2D a t a . An important property of the fractional 
Fokker-Planck equation is that in the long time limit 
Boltzmann equilibrium is obtained [32j 



p cq (x) = 



exp 



V(x 
T 



z 



(38) 



provided that the force is binding. Recently numerical 
methods which give the sample paths of the fractional 
Fokker-Planck equation were investigated in detail [2(| 
H3, [H, [H, H3| • Such paths or the corresponding CTRW 
trajectories yield non ergodic behaviors [lj| . For example 
in [16[ the Lamperti PDF Eq. ([23]) of the occupation 
fraction was obtained from the fractional equation (|37[) . 
However so far distributions of time averages of physical 
observables were not considered in detail. 

We investigate the time average of the observable O = 
x with — oo < x < oo so we are dealing with a continuum 
situation. Taking the continuum limit of Eq. (TH 
find 



1 J °° 

f a (x) = lim Im^£ 

7T e^Q J_ 



dxp eq (x) (x — x + ie) a 
' dxp c<1 (x) (x — x + ie) a 



(39) 



which for < a < 1 is rewritten as 



fa 0) = 



sin na 



/<_ 1 (x)7>(S) + />_ 1 (5)7<(x) 



m 2 



2 cos iral> (x) I< (x) 



(40) 



and similarly for {x) and (%)• When a — > 
we have lim Q ^o fa {%) = p cq (x) which is the continuum 
limit of Eq. (|21|) . In the ergodic limit a — > 1 we find 

h(x)=S(x-(x)). 



A. Free particle. 

As an example consider a particle in a domain —1/2 < 
x < 1/2 undergoing an unbiased fractional random walk 
with reflecting walls. This is a free particle in the sense 
that no external field is acting on it. The time average 
of the particle's position x is considered, and obviously 
for this case p eq (x) = 1/1 for —1/2 < x < 1/2. Using Eq. 
(HI we find the PDF of x 



fa (x) 



N, 



a (4 



T- 



l 



1 2(l+oi) 



+ 



,2(l+a) 



(!)" 



l+a 



cos ira 



( 43 ) 

where N a = (1 + a) sin7ra/(7ra). When a — > 1 we have 
ergodic behavior f a (x) = S(x) since (x) = while in the 
opposite limit f a ^o (x) = 1/1 for \x\ < 1/2 which is the 
uniform distribution, reflecting the mentioned localiza- 
tion of the particle in space when a — > 0. 



Harmonic Oscillator 



We consider the time average x of a particle in a 
Harmonic force field V(x)/T = x 2 so that p cq (x) = 
exp(— x 2 )/y / 7T and (x) = 0. Using Mathematica the 
integrals Eqs. (|41l42p can be calculated explicitly and 
expressed in terms of tabulated confluent Hypergeomet- 
ric functions. In Fig. [T] the PDF of x Eq. (J40J) is 
presented, and a transition from a narrow distribution 
when a — > 1 to a Gaussian distribution when a — » 0, 



lim Q ^o fa {%) = P° q (x) is found, 
is easy to show that (x) = 0, 



Using Eq. (gDJ it 
(x 2 ) = (1 - a)(x 2 ) 



with (x ) = 1/2 for p cq {x) under consideration, and 
(x 4 ) = (1 - a) (3 - 2a){x 2 ) 2 . Only when a -> we 
have Gaussian statistics with (or) = 3(x 2 } 2 . The PDF 
faix) at its maximum on x = is 



fa (X = 0) 



r(f)tan(^) 



r( 



l+a-. 
2 / 



(44) 



FIG. 1: The PDF of a; for a particle in a Harmonic force field 
Eq. (|40p . We find a transition between an ergodic behavior: 
a delta distribution of a; when a — > 1, to the localization limit 
where the distribution of x is Gaussian when a —* 0. 

which is equal to p eq (x = 0) = 1/\/tt when a — > and 
diverges when a — > 1 as expected from an ergodic behav- 
ior. For the Harmonic oscillator and the Free particles 
the maximum of f a (x) is found on the ensemble average 
(x) = 0, so the most likely result for x is (x). In the next 
subsection we consider a case where a minimum of f a (x) 
is found on the ensemble average. 



C. Double well potential. 

An interesting case is the symmetric double well po- 
tential V(x)/T = (x 4 /4 - x 2 /2)/T so (x) = 0. When 
T — > 0, p eq (x) has two peaks centered on the two local 
minima of the double well potential. In this low temper- 
ature case and in the limit ct-*0we expect to find the 
particle either in the left well or in the right well for a 
time scale comparable to the measurement time. Hence 
when T — > and a — > the PDF of the time average x, 
f a (a?) is a sum of two delta functions since either x — 1 
or x = — 1. When a — > 1 we expect an ergodic behav- 
ior, and then PDF /i (x) = 6(x), since (x) — 0. So for 
low temperatures we expect a transition in the behavior 
of f a (x) from a bimodal shape when a — ► to a PDF 
with a single peak centered on zero when a — ► 1. Hence 
we will have a critical value a c . For cv < cv c the shape 
of / Q (x) is bimodal with a minimum on x = 0, while 
for a > a c a maximum on x = is found. These low 
temperature behaviors are shown in Fig. [2] For high 
temperatures (compared with the barrier height) the bi- 
modal solution of p eq (x) turns into a flatter shape. Since 
lima^o fa {%) — p cq {x) we will not observe the bimodal 
shape of lim^^o f a (x) when T — + oo. Such high temper- 
ature behavior is shown in Fig. [3] 

Investigating the extremum of p cq (x) on x — it is 
easy to show that a c is finite for any finite temperature 



FIG. 2: The PDF of x for a particle in the double well 
potential with T — 0.01. A transition between a bimodal 
behavior when a < a c to a PDF with a peak on x = when 
a > ct c is observed (a c — 0.59 for this case). In the ergodic 
limit a — > 1 f a (a?) is a delta function centered on the ensemble 
average (a;) = 0. In the localization limit a — > f a (x) is equal 
to the population density p cq (aT). 

and a c — > when T — > oo. For T — we have only 
two states in the system, either x — ~\ or x — 1 cor- 
responding to two minima of the double well potential. 
The analysis is then very similar to the two state ballistic 
Levy walk model [HI [38[ . Clearly x is the residence time 
in state x = 1 minus the residence time in state x = — 1 
divided by the measurement time t. Since the sum of 
these two residence times is the measurement time we 
can use the Lamperti distribution Eq. (|23[) to predict 
the distribution of x. So when T — > 

lim f a (x) = 

2sin7ra (l-x 2 ) a 1 

TT (1 + x) 2a + (1 - x) 2a + 2 (1 - a; 2 )" cosTra 

(45) 

which was found already in [2(| [3t|. We find that 
a c = 0.59461 ■ • • when T — > 0. The behavior of cv c versus 
temperature is shown in Fig. [4] and the transition be- 
tween the low and high temperature cases is presented. 

D. Possible Physical Applications 

It is interesting to verify in experiments our theoretical 
predictions and here we discuss three examples. Gener- 
ally systems with CTRW type of dynamics are natural 
candidates for the investigation of weak ergodicity break- 
ing, provided that information of single particle dynamics 
can be recorded. 

Sub- diffusion {x 2 } ~ t a of a bead in a polymer network 
was measured by Wong et al [4(j. The measured [4(| 
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X 



FIG. 3: The same as Fig. © however now T = 7. The 
bimodal shape presented in Fig. ([2]) is smoothed and we 
barely observe bi-modal behavior, since the equilibrium den- 
sity p eq (x) is not centered around the two minima of the dou- 
ble well potential, when the temperature is high. 
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FIG. 4: The critical exponent a c versus temperature for a 
particle in a double well potential. a c marks the transition 
between a local minimum to a local maximum of the PDF of 
x on x — 0. When T — > oo a c — > since in this limit and 
when a —> f a (x) is equal to the population density p eq (x) 
which does not "feel" the wells when the temperature is high 
(single peak). When T^Owe have a c ~ 0.595 . 



exponent a depends on the ratio of the size of the bead 
and the linear size of the mesh of the network I (roughly 
a /im). We suggest to add an external binding field, for 
example an harmonic trap. The time averages of a single 
particle coordinate can then be measured, and accord- 
ing to our theory its distribution is given by Eq. (|13[) . 
Such measurement can provide insight into the nature of 
disorder, for example is it quenched or annealed [22| . 

Messenger RNA molecules inside live E. coli cells ex- 
hibit anomalous diffusion (a; 2 ) ~ t a and a = 3/4 fill ]. 
Due to the finite size of the cell the motion is bounded. 
It would be interesting to investigate time averages of 
the position of the single molecule, or occupation time 
statistics, to investigate deviations from ergodicity. Our 
theory gives a prediction for the distribution of these ob- 
servables, which can be tested in experiment. 

Blinking quantum dots exhibit ergodicit y b reaking 
which is already measured in experiments [^.llil. H3|. So 
far a simple two state picture of the quantum dots was 
used, either the dot is on and it emits many photon, or it 
is off [lfj . Then ergodicity breaking of the time averaged 
fluorescence intensity is similar to the time average of a 
particle in the double well potential in the limit of T — > 0, 
in the latter case either the particle is the left well or in 
the right well. More careful analysis reveals that some 
dots deviate from a simple two state process Then 
our more general theory can be used in principle to pre- 
dict distribution of time averaged intensity beyond the 
existing two state approach. 



V. FROM CONTINUOUS TIME RANDOM 
WALK TO WEAK ERGODICITY BREAKING 

In this Sec. we derive our main results using the 
CTRW approach. To reach Eq. ([T5)) we assumed among 
other things that: (i) the PDF of the occupation time 
ti is the one sided Levy PDF Eq. ^ and (ii) that 

the total measurement time t — 53j=i b ^ s a random 
variable, while in Physical experiments the measurement 
time is fixed. These assumptions are relaxed now using 
two types of CTRW models. Thermal CTRW describe 
a Physical situation where the particle is undergoing the 
random walk in contact with a thermal heat bath. In 
this case the equilibrium distribution of an ensemble of 
particles is Boltzmann's distribution. The second case 
describes a system far from thermal equilibrium, where 
a non-thermal equilibrium is reached. 

We consider a renewal process for a system with L 
states i = 1, L. The system starts in state i, it waits in 
this state until time t±, it then jumps to some other state 
say state /, it waits in state I until time ii and then makes 
another jump. The sojourn times between jump events r 
are independent identically distributed random variables 
with a common PDF ip( T )- Our focus is on the case 
where ip{ T ) has a long tail ip{ T ) ^ t~( 1+q ) when t — > oo 
so (r) = oo when < a < 1. After waiting in a state i 
a transition to state j ^ i takes place, with probability 

Wji (0 < Wji < 1 , J2j=i w ji = 1 and w u = 0)- We 
assume that transition probabilities Wji are such that in 
the limit of long measurement times all states are visited 
what ever is the initial condition. In other words the 
system is not decomposable into non-accessible regions 
in the space it samples, where once the system starts in 
a certain region it cannot explore all other states. Such 
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a case corresponds to strong ergodicity breaking. 

Let N be the random number of jump events (re- 
newals) in the time interval (0,i). Dots on the time axis 
on which jumps from one state to another happen are 
denoted with and clearly < t < tjv+i- Let n% 
be the number of transitions out of state i and clearly 



Let Tj be the I th sojourn time in state 



i. And let k be the state of the system at time t. A 
schematic presentation of the process with three states 
is shown in Fig. [5l Statistics of the number of renewals 
N in (0, t) is a well investigated problem [H, [H, |42j for 
example (A) ~ t a . 

For the CTRW on a one dimensional lattice the states 
i = 1, ...L correspond to the position of the particle on a 
finite lattice. Then for example the time average of the 
coordinate of the particle is x = where ti is 

the occupation time in state i (see Fig. [5]). However our 
considerations are mo re g eneral. For example for blinking 
quantum dots [H, [n| El, E3] one state (say i = 1) may 
denote an on state in which many photons are emitted 
and state i = 2 is the off state. This system is non- 
thermal since it is driven by a strong laser field. On 
the other hand the CTRW dynamics of a probe bead 
immersed in a polymer actin network [40l | is an example 
for a thermal CTRW motion in a system with a well 
defined temperature T. 

A specific example is the CTRW on a lattice with 
jumps to nearest neighbors only. A particle on i has a 
probability of jumping left qi and a probability of jump- 
ing right 1 — qi so 



-l,« 



and 



= 1 



(46) 



Reflecting boundary conditions qi, = 1 and q\ — are 
assumed. For i ^ 1, L we must have qi ^ and qi ^ 1 so 
that all lattice points be visited. Between jumps the par- 
ticle waits on a lattice point. The waiting times between 
the jumps are independent, identically distributed ran- 
dom variables with a common PDF ip(r). This type of 
random walk leads to anomalous subdiffusion (x 2 ) ~ t a 
when qi = 1/2 and the system is infinite 0, 52 . 

The ratio i?j = rii/N is called the visitation fraction. 
The population fraction p eq is found by considering the 
ensemble of M non interacting systems. Letting these 
systems evolve from some initial condition and waiting 
for the long time limit, limA/^oo Mj/M = p° q where Mi is 
the number of systems in state i. The population fraction 
is determined from w ■ p cq = p eq . 

After many jumps N — ► oo and for any initial condition 
the visitation fraction reaches an equilibrium and 



P 



eq _ 



= lira 

N-xx 



(47) 



so w ■ v — v. To see this note that the visitation fraction 
is given by = J2 n =i ®i( n )/N '•> where n is a counter of 
the number of jumps, and 6i{n) = 1 if the particle is 
on i after n steps, otherwise it is zero. In the long time 
limit the PDF of Vi will converge to a narrow distribution 
centered around its mean so Vi ~ En=i(^i( n ))/-^- Let 



Pi(n) be the probability to be on i after n steps. By 
definition Qi (n) = 1 with probability pi (n) and 6i (n) = 
with probability 1 — Pi (n) . Hence 



E^=iP<( n ) 



N 



(48) 



(vi,-~,Vl), P{n) 



or m vector notation v 

(• • • ,Pi(n) • • •) we have v ~ Sn=i P( n ) I 'N ■ This means 
that ergodicity holds in discrete time, where the oper- 
ational time is the number of steps, not the real time. 
Hence the term weak ergodicity breaking Q is very ap- 
pealing. Multiplying Eq. (|48|) with w from the left and 
using p(n+l) — wp(n) we have w-v ~ J2n=i 
Hence when N — > oo we have w-v = v which holds in the 
long time limit. It is important to realize that the visita- 
tion fraction and the population fraction are equal since 
all sojourn times have a common distribution 4>(t). Wc 
will later consider the more general case where different 
states may have different waiting times PDFs. 

For the one dimensional CTRW on a lattice the equi- 
librium population and hence the visitation fraction is 
determined from Eq. (|46[) 



pr = ^iP^ + (l-9i-i)P^i. (49) 
Using reflecting boundary conditions and Eq. (|49[) 



P 



eq 



lim Vi 



1 



AT— »oo 

and from normalization 

L 



-ni. 



1 - qi qk 



1 - <lk cq 

- "Pi 



eq eq 

Pi = V 



l= 2 1 qi 



(50) 



(51) 



When the particle undergoing the CTRW process is cou- 
pled to a thermal heat bath, we apply usual detailed bal- 
ance condition on the transition probabilities [l5| . In this 
case the visitation fraction will be described by Boltz- 
mann statistics. For example if qi is a constant q > 1/2 
the random walk is biased, which Physically corresponds 
to an external force field F < driving the particles to 
the left. Using lattice spacing a and letting the system be 
semi- infinite L — > oo, thermal detailed balance condition 
gives the ratio between the probability of jumping left 
from point i and the probability of jumping right from 
point i — 1 



i = i = ex P (hfr- ■ 

1 - 1 - q T 



(52) 



Using Eqs. (|47I50I51I52|) Boltzmann's statistics holds 
both for the visitation and the population fractions 



exp (-jf) 
Z 



(53) 



for i > 1 where Ei = \F\ai is the potential energy, Z is 
a normalization and for the reflecting boundary p eq = 



10 



1 = 1 l = 2 l = 3 



FIG. 5: A schematic diagram of the process for a system 
with three states, starting in state 2 and ending in state k = 
1. In this example the occupation times are ti = t — tz, 
t 2 = (ti - 0) + (h - i 2 ) = n 2 + r| and t 3 = t 2 - ti = rf. 



[1 — exp(— |F|q/T)1 /2. More general thermal detailed 
balance conditions [15J show that Eq. is valid for 

binding force fields and not limited to the case F being 
a constant. 

The time average of a physical observable is as before 



q _ J2i=i Oik 



(54) 



where now the measurement time t is fixed and Ej— i ti 



t. Let O = {Oi, • • • 0_l}. We consider the moment gen- 
erating function, 



ft,d ( u ) = ( ex P ^ 
and in double Laplace space 



(55) 



e s *(exp 



)dt 



(56) 



so s, £ and it, Ei=i O^j are two Laplace pairs. Let 
n = {ni, • • • , til}. We consider the generating function 
conditioned that the system made N transitions and ft 
describes the number of renewals in each state. The oc- 
cupation time in state i ^= k is 



i=i 



(57) 



and for state k 



t - t 



N 



ran, 

\ ~" k 



1=1 



(58) 



The time t — t jy is called the backward recurrence time, it 
is the time between the last jump event in (0, t) and the 
measurement time t. Using Eqs. (|57I58|) the conditioned 
generating function is 



fs.d.N.n (") = ( / d * eX P 



-St-U0 k [t-t N + J^T? ) - U °iJ^ T l 

i=l^k 1=1 



1=1 



I (tjy < t < t N+ i)) (59) 



where I(x) = 1 if the condition in the parenthesis is true, other wise I(x) — 0. First we integrate over t and obtain 

-sijv p — sijv + i— uOk(tN+i — tjv) 



fs.O.N.n ( U ) (" 



s + uOu 



(60) 



then we use the assumption of independent and identically distributed sojourn times r, and the identities tn 
£f=i Er=i ^. ftv+i = t N + r^ +1 to find 



nf =1 ^ ni (s + i-v(* + «0k) 



(61) 



r 



where ip(s) = J °° ip{ T ) ex P( — st)At is the Laplace trans- the system will reach an equilibrium for the number of 
form of ip{ T )- 

In the limit of long measurement time, corresponding 
to the usual small s and u limit, their ratio being finite, 



renewals in each state. Namely from Eq. (|47|) the visita- 
tion fraction will satisfy 



Vi = lim — = p° 

TV— oo N 



(62) 
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FIG. 6: Trajectory of a single CTRW particle on a lattice, 
with q = 0.3 (solid curve, blue online). Long waiting times, of 
the order of the measurement time, dominate the landscape. 
The time average a; is a random variable (dotted red curve). 



We use Eq. (162I) in Eq. ([61]) . then insert the usual small 
s behavior [aEl [H, El 

tp(s) ~ 1 - As a (63) 
which corresponds to ip( T ) ~ ^4r _ ( 1+Q ^/|r(— a)\ and find 



(64) 



Summing over all final states k, with the final weights, 
we find 



fs.o ( u ) 



(65) 



Using a method found in the Appendix of Ref. [26| we 
invert Eq. (|65p and find our main Eq. ([Tc 



FIG. 7: Same type of non-ergodic behavior as found in Fig 
[5] however now a = 0.75. 




x 10 



FIG. 8: Same as Fig. [6] however now ip(r) — 5t~ 6 for r > 1 
so a = 5. After a short transient the time average converges 
to the ensemble average on 15.5 indicating ergodicity. 



A. Numerical Examples 

To demonstrate our results we consider an unbiased 
CTRW. We consider a model with L = 30 sites on i — 
1, ...30 jumps are to nearest neighbors only with qi = 1/2 
with periodic boundary conditions. We used the waiting 
time PDF ip(r) = ar^ ^ 1 for r > 1 otherwize tp( T ) = 0. without fitting. In Eq. (fT3|) we use p° q = 1/L which is 



found for strong ergodicity breaking. For a waiting time 
PDF with a = 5, we have a finite average waiting time, 
and hence as shown in Fig. [5]we find an ergodic behavior 
x ~ (x) = (L + l)/2 = 15.5 in the long time limit. 

In Figs. [gnUl we present the PDF of x for a = 0.75 
and a — 0.3 respectively. Comparison between simu- 
lations and theory Eq. (fTB")) show excellent agreement 



Simulating trajectories of a single particle we calculate 
the time average x and then repeat the experiment many 
times and construct the distribution of x. 

In Figs. 16171 we show single particle trajectories and 
their time average for a — 0.3 and a — 0.75 respectively. 
The time average is clearly a fluctuating random variable, 
due to long sticking times in states of the system. Notice 
that the particle visits all lattice points, and the phase 
space is not decomposable into inaccessible regions as 



the obvious population probability. The number of real- 
izations was 120000 and the simulation time t = 10 8 , 10 12 
for Figs. EHO] respectively. In Figs. I9|10l we also show the 
continuum approximation Eq. (|43p . From Figs. I9|10l we 
see that the structure of the lattice is encoded in the dis- 
tribution of the time average x. Since the observable at 
any given moment of time attains the values x — 1, ...30 
we have 30 divergences in Figs. 191101 in agreement with 
the more general rule discussed under Eq. (|16p . When 
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FIG. 9: Distribution of the time average x for a = 3/4. Nu- 
merical simulations of the CTRW process on a lattice (dots) 
are well approximated with the continuum limit Eq. (|43[) 
(solid curve). The dotted curve with divergences on the lat- 
tice points is the analytical theory Eq. (|13|) . 
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FIG. 10: Same as Fig. |9]for a = 0.3. Now the fluctuations 
are larger compared with the a = 3/4 case and the underlying 
structure of the lattice is important. Ofcourse bin size must 
be made small enough for the lattice effect to be observed. 



the system is made large these effects become negligible 
and the continuum approximation works well. 



B. Non identical waiting time distributions 

The CTRW considers a situation where a single waiting 
time PDF ip(r) describes the dynamics. What happens 
when the waiting times in the states i are not identi- 
cally distributed? For ex amp le consider the mentioned 
blinking quantum dots [§, O EH E2 ■ The quantum dot 
when interacting with a laser field will switch between 
an on state (+) where many photons are emitted and an 



off (— ) state. The sojourn times in state on (and off) 
are independent identically distributed random variables 
[43l ]. The PDF of sojourn times in states on and off fol- 
low power law statistics atleast within a long experiment 
time [l2j and in Laplace space ip±(s) ~ 1 — A±s a + ■■ ■. 
This two state renewal process is characterized with two 
amplitudes A+ and A- and in this sense it differs from 
the usual CTRW. It is worthwhile noting that the visita- 
tion fraction in states ± clearly satisfy 



lim 

N^oc 



v+ = lim — = 

N-*oo N 



n_ 1 
lim — = — 

N^oc N 2 



(66) 



in the limit of long measurement time, where N is the 
total number of transitions between states on and off. If 
we consider an ensemble of M independent blinking dots 
the population fraction of the number of dots in state on 
(+) M+ and off (-) M_ is 



p+ = lim — — - 



A± 



A_ 



A, 



(67) 



where the population fraction is measured in the limit of 
long measurement time. So here the visitation fraction 
and the population fraction are non-identical if A + 7^ 
A_. 

More generally we consider the renewal dynamics with 
power law waiting times in each state however now 



tPi(s) ~ I — Ai 



(68) 



when the Laplace variable s — > 0, Ai > for i = 1, • • • , L 
and < a < 1. In this case the population fractions are 
related to the visitation fractions according to 



lim 



N- 



(69) 



and w ■ v = v. Now the main Eqs. derived in this Sec. 
must be modified, for example Eq. (|6ip 



f s.O.N. n ( u ) 



(70) 

Then using Eqs. (|68|69|70p we find Eq. (JT3J) (the method 
is nearly identical to the case where all the waiting times 
are identical). Thus while the dynamics clearly differs 
from the usual CTRW with a single waiting time PDF, 
and the visitation fraction is not identical to the popu- 
lation fraction, our main Eq. for the distribution of the 
time averages Eq. flT5|) is still valid. 

Another situation is when the system has different 
types of waiting times, for example some states may have 
exponential waiting times while others follow power law 
statistics. Or we may have some states with a power law 
waiting time PDF with an exponent < a\ < 1 and 
for other states an exponent < a<x < oti. Also in this 
case our main result Eq. (flU)) will be valid. In the long 
time limit the system will occupy only the states with 
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the smallest a < 1. Only those states are relevant for 
the calculation of the distribution of time averages Eq. 
(fT"3|) . Other states might be visited many times so their 
visitation fraction is not necessarily small, still the time 
the system spends in these states is short and they do 
not contribute to the time average. 

VI. DISCUSSION 

To summarize we have obtained very general distribu- 
tions of time averages of physical observables of weakly 
non-ergodic systems Eqs. (|13I39[) . Unlike usual ergodic 
statistical mechanics where the time averages are equal 
to the ensemble averages, we find large fluctuations of 
time averages. Unlike strong ergodicity breaking, the 
space is not separated into inaccessible regions. Instead 
the system explores all states and the number of visits 
per state i rn is large. However due to the power law 
sojourn times the dynamics is non-stationary and non- 
ergodic. Since exploration of space is possible in weak 
ergodicity breaking, we were able to construct a rather 
general statistical theory for the distribution of time av- 
erages which is valid in the long time limit and does not 
depend on the initial conditions of the system. Further, 
the exploration of the cells i — 1, ...L leads to relations 
between the non-ergodic statistical properties of a single 



system, and the equilibrium populations of an ensemble 
of systems. Hence the distribution of time averages Eq. 
(fl"3f depends on population probabilities p° q s and when 
thermal detailed balance is applied these probabilities are 
given by Boltzmann's canonical law. Such behavior was 
found in two classes of models: (a) for systems with a 
common power law waiting time PDF ip{ T ) where the 
visitation fraction is equal to the population fraction and 
(b) for systems with different waiting times PDFs where 
the visitation fraction is not equal to the population frac- 
tion. In both cases the population fraction is not iden- 
tical to the occupation fraction, the latter remaining a 
random variable when < a < 1. Due to the large 
number of applications of the CTRW model, our theory 
may find its applications in several systems. The math- 
ematical foundation of the theory is Lamperti's general- 
ized arcsine law and Levy's statistics, while usual type 
of statistical mechanics is based on the Gaussian central 
limit theorem. Due to the deep relations between the 
stochastic CTRW model and other models of anomalous 
diffusion, e.g. the quenched trap model and determin- 
istic dynamics our non-ergodic theory might find more 
general justification. 
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